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Abstract 

It was found by Perdew, Parr, Levy, and Balduz [Phys. Rev. Lett. 49, 1691 (1982)] and by 
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Sham and Schluter [Phys. Rev. Lett. 51, 1884 (1983)] that the exact Kohn-Sham exchange- 
correlation potential of an open system may jump discontinuosly as the particle number crosses 
an integer, with important physical consequences. Recently, Sagvolden and Perdew [Phys. Rev. 
A 77, 012517 (2008)] have analyzed the discontinuity of the exchange-correlation potential as 
i the particle number crosses one, with an illustration that uses a model density for the H~ ion. 

In this work, we extend their analysis to the case in which the external potential is the simple 
harmonic confinement, choosing spring-constant values for which the two-electron hamiltonian 
has an analytic solution. This way, we can obtain the exact, analytic, exchange and correlation 
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the particle number crosses one without introducing any model or approximation. We also discuss 
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potentials for particle number fluctuating between zero and two, illustrating the discontinuity as 
the particle number crosses one withe 
exchange and correlation separately. 
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I. INTRODUCTION 



Kohn-Sham (KS) density functional theory (DFT) (see, e.g., pQ) is a successful method 
for electronic structure calculations, thanks to its unique combination of low computational 
cost and reasonable accuracy. In the Kohn-Sham formalism, the total energy of a many- 
electron system in the external potential V ext = J2i v ext( r i) is rewritten as a functional of 
the one-electron density p(r), 

E[p] = T s [p] + J drv cxt (r) p(r) + U[p] + E xc [p}. (1) 

In Eq. (JTJ, T s [p] is the kinetic energy of a non-interacting system of fermions (usually called 
KS system) having the same one-electron density p of the physical, interacting, system. The 
Hartree energy U[p] is the classical repulsion energy, U[p] = | j dr j dr'p(r)p(r')\r — r'| _1 , 
and the exchange-correlation functional E xc [p] must be approximated. Minimization of 
Eq. (jlj with respect to the orbitals forming the KS determinant lead to the KS equations. 
Thus, instead of the physical problem, in KS DFT we solve the hamiltonian of a model 
system of non- interacting fermions in the one-body local potential Vks = ^2i v Ks( r i), with 

fKs(r) = fext(r) + vn(r) + v xc (r) (2) 

VH(r) = / dr'^ (3) 
5 EM 

Vxc[ } Sp(r) ' 1 ] 

and we recover the energy of the physical, interacting system, through the sum of the two 
functional U[p] + E xc [p\. 

In Refs. [21 El H] an analysis of KS theory for systems with fluctuating particle number 
lead to the conclusion that the exact exchange-correlation potential f xc ( r ) may jump dis- 
continuously by a spatial-independent constant as the particle number crosses an integer, 
with important physical consequences. In the last years there has been new interest in the 
derivative discontinuity of E xc [p] (see, e.g., El [3, EJ |9|, UHl [TT]), and its existence has been 
questioned in Ref. [T2] . In a recent paper, Sagvolden and Perdew [13j have given further 
support to the assumptions used to find the exchange-correlation potential discontinuity, 
and they have illustrated the discontinuity when the particle number crosses one, using a 
model density for the H~ ion. They have also rigorously proved that the von Weizsacker 
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functional, 



TMP] = \j dr\V^)\ 2 (5) 



is the correct T s [p] for a system with particle number N < 2. 

In this work we repeat a similar analysis for electronic systems with particle number 



fluctuating between zero and two when the external potential is harmonic, u ex t(r) = \uj 2 r 2 . 
Taut [H] has shown that the corresponding hamiltonian for N = 2 electrons is analytically 
soluble for some special values of u, which means that, in such special cases, we can calculate 
analytically the exact interacting density and the exact exchange-correlation potential [15J. 
This way, we can illustrate the derivative discontinuity without relying on any approxima- 
tion. The paper is organized as follows. In Sec. [IT] we report the equations used to extract 
the exact exchange and correlation potentials from Taut's [H] analytical solutions. The 



corresponding results are shown and analyzed in Sec. |III[ and the last Sec. [IV] is devoted to 
concluding remarks. 

II. THEORY 

We consider a system with particle number fluctuating between and 2, in the harmonic 
external potential f ex t( r ) = \oj 2 t 2 . This system can be thought as arising from identical 
quantum traps whose centers are separated by a very large distance [16J. If the energy 
of the physical system is convex as a function of the integer particle number M, i.e., if 
E(M) < (E(M + 1) + E(M — l))/2, then the system with non-integer particle number N 
is an ensemble of only the two systems with integer-particle number M and M + 1 such 
that M < N < M + 1 [2]. In the case treated here, t> cx t( r ) = \uj 2 r 2 , the energy of the 
non-interacting system is convex, and the electron-electron interaction seems to make it 
strictly convex (see, e.g., the results for two-dimensional harmonic traps of Refs. [T7] or the 
three-dimensional case treated in Ref. [18J). The density of the system with particle number 
< N < 2, p N (r) = p N (r), is then equal to E3 

, , , N Pl (v) 0<N<1 
PN{r) = { (6) 
(2 - iV)pi(r) + (N — l)p 2 (r) 1< N < 2 
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The exact Kohn-Sham potential is, up to a constant, given by the functional derivative with 
respect to Pat(t) of the von Weizsacker functional of Eq. (J5|, 

VKsir) = — + const. (7) 

The Hartree potential vu(r) can be easily calculated by plugging Eq. ^ into Eq. 
The exchange and correlation potential v xc (r) is then calculated from Eq. Q, v xc (r) = 
v Ks( r ) — — f ex t( r )- Setting v xc (r — > oo) = for any fixed N determines the arbitrary 
constant in the exchange-correlation potential. 

Splitting the potential v xc (r) into exchange and correlation is subtle. When < N < 1 
we know that v x (r) must exactly cancel the Hartree potential, v x (r) = —Vu(r), and that, 
since the correlation energy is zero, v c (r) = 0. For 1 < N < 2 the sum U[pn] + E x [p n ] should 
be equal to the expectation of the electron-electron repulsion V ee over the non-interacting 
ensemble density matrix that yields Pn{^)- The non- interacting density matrix that yields 
PAr(r) and corresponds to the kinetic energy functional of Eq. ^ is 



T = (2 - A0|$i)($i| +(N- 1)|$ 2 )($ 2 |, (8) 



with 



•.M = JfHp (9) 



Thus 



•»<™) = JmpL^EnpiL. ( 10 ) 

U\p N ] + E x [p K ] = Tr(r„l>«) = 2( - N N ~ l ' u[p N ]. (11) 

Unlike the Hartree functional U[pn] and the non-interacting kinetic energy functional 
T v w[pn] of Eq. ([5]), which have a simple explicit dependence on P7v(r) alone, we see from 



Eq. (11) that the exchange functional E x \p^\ also explicitly depends on the particle number 



N . If we take the functional derivative 5/5pN(r) of both sides of Eq. (11) at fixed constant 
N, the exchange potential is not discontinuous at N — 1. If, instead, we allow little varia- 
tions of iV with respect to pjv(r), by taking into account that N = J drp^(r), we obtain, 
for the whole range < iV < 2, 

, W = l "" H(r) ° " N ~ 1 (12) 

xK ' ' (N 2 -2N+2) i \ . 2(2- N)--- - 



^ va(r) + £ W ± U[ PN ] 1<N<2 



Equation (12) shows that, within the definition of exchange of Eq. (11) and allowing vari- 
ations of N with respect to p, the discontinuity of the exchange potential at iV = 1 is a 
spatially-independent constant equal to 

v x (r)\ N ^i+ ~ v x {r)\ N ^i- = 2U\px\. (13) 

Exact exchange in an open system of fluctuating electron number as also been widely dis- 
cussed in Ref. [19], where only the energy, and not the potential, has been analyzed. 

When the external potential is harmonic, v ext (r) = ^u 2 r 2 , we have pi(r) = 
cg> 3 / 2 /7t 3 / 2 e~ wr2 . The ground-state wavefunction of the hamiltonian with N = 2 electrons 
has the form 

#(ri,r 2 )=77(i2)V>(ri 2 ). R= ^ * r 12 = |r 2 - n|. (14) 

The center-of-mass wavefunction r)(R) is a simple ground-state three-dimensional harmonic 
oscillator state. Taut [T4] has shown that for some special values of u also the relative 
wavefunction V( r i2) has an analytic form. For some of these special u- values we have 
calculated the corresponding electronic density p 2 (r), 

_ — LOT 2 POO 

p 2 (r) = C- n / dr l2 r 12 t(r 12 ) 2 ( e -(— ) 2 - e ^)') , (15) 

2ur J V / 

1 2 

where t(x) = e^ x ip{x), and C is a normalization constant. The function t(x) is a polyno- 
mial of some finite order n (depending on u), so that the corresponding p 2 (r) is completely 
analytical (see Appendix [AJ . This allows us to obtain the exact exchange-correlation poten- 
tials up to any large r exactly, without introducing any approximation. The large-r part of 
the potentials is crucial in order to illustrate the derivative discontinuity. 

Notice that similar calculations of the exact exchange-correlation potential and of other 
DFT-related properties for the N = 2 system have been carried out in the past, mainly 
aimed at comparing the exact potentials with current approximations (see e.g., [T5l [20| [2Tj ) . 
Here our aim is to calculate the exact potentials for < iV < 2 in order to study their 
discontinuity as the particle number crosses one. 

III. RESULTS 



We considered three cases for which ^(r^) of Eq. (14) is analytical: to = 1/2, u = 1/10 



and to = (35 — 2v57)/712 w 0.0173 (see Appendix |A[). The exact exchange-correlation 
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FIG. 1: The exact exchange-correlation potential for a system with particle number N fluctuating 
between zero and two in the external potential ^ojr 2 , for three values of u. The magnitude of 
the discontinuity at N = 1 is equal to E2 — 2E\. We have for uj = | E2 — 2E\ = |, for 00 = 
E 2 -2E l = \, and for u » 0.0173 E 2 - 2E 1 = 4w w 0.06938. 
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potentials for N slightly below and slightly above 1 are reported in Fig. [T} The qualitative 
behavior is similar to the one reported by Sagvolden and Perdew [T3] obtained from the 
model density of the H~ ion: as iV — > 1 + the exchange-correlation potential is, over a 
larger and larger range, more and more equal to the N — > 1~ potential plus a constant. 
The magnitude of the constant is, as expected [21 El H], I — A, where I is the ionization 
potential and A the electron affinity For closed-shell two-electron systems this amounts to 
I — A = E 2 — 2Ei, where E 2 is the energy of the system with N = 2 electrons and E± is the 
energy for N = 1. Figure [I] clearly shows that, as argued in Refs. [21 El IH H3], for any fixed 
given N (above or below 1) 

lim v xc (r) = 0, (16) 

but 

lim lim v xc (r) = I — A. (17) 

r— >oo JV— »1+ 

Notice that when u = 1/2 and u = 1/10 the N = 1.0000000001 v xc {r) starts to decay to zero 
for distances much larger than those considered in the figure. The case of the ensemble H and 
H studied by Sagvolden and Perdew [13] resembles most closely to the last case considered 
here, uj ~ 0.0173. This is due to the fact that the N = 2 system is considerably 
correlated. In the case treated here, f cx t(^) = \oj 2 r 2 1 the N = 2 system becomes more and 
more correlated as uj — > 0. This is also consistent with our previous work [22], in which we 
compared the dependence on iV (in the range 1 < iV < 2) of the non-interacting kinetic 
energy of the He series with the one of the Hooke's atom series, finding a resemblence of the 
H~ case (He series) with the low-u; case (Hooke's series). 

IV. CONCLUSIONS 

The discontinuity of the exact Kohn-Sham exchange-correlation potential was argued in 
Refs. [21 El H] , and is known to have very important physical consequences. In particular, the 
absence or the underestimation of this discontinuity in local and semilocal functionals ex- 
plains why they often produce a qualitatively incorrect dissociation limit for non-equilibrium 
nuclear positions and why they underestimate charge-transfer excitation energies in time- 
dependent DFT. In this work we have calculated the exact Kohn-Sham exchange and corre- 
lation potentials as the particle number crosses 1 for a system which is analytically soluble. 
This way, we have illustrated, at least in one case, the existence of this discontinuity without 
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relying on any approximation, thus providing further support to the assumptions that were 
used for its prediction. 
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APPENDIX A: ANALYTICAL DENSITIES FOR N = 2 



As shown by Taut [H], for a set of special values of u the ground-state wavefunction of 
the two-electron hamiltonian has the form 



^(n,r 2 ,r 12 ) = Ce~^ r * +r *h{r 12 ), 



(Al) 



where r\ = |ri|, r 2 = |r 2 |, r\ 2 = |r 2 — ri|, t(x) is a polynomial of finite order n (depending 
on u), and C is the normalization constant. The easiest way to obtain the density is to use 
the coordinates of Coulson and Nielson [23], r±, r 2 and r 12 , for the integration over r 2 , 



-u)r roo 



P2{r)=C- 



dr u r 12 t(r 12 ) 2 I dr 2 r 2 e~ 

<> J\r2-ri2\ 



\r2+ri2\ 



where C is again a normalization constant. Integration over r 2 yields Eq. (15) 



(A2) 



1. w = l/2 

The density for uj = 1/2 has been reported by several authors (see, e.g., Ref. [21]), 

4 (r 2 + 1) erf 



e 2 



4 (8 + 5^r) 7T 3 / 2 



P2{r) 



where erf(a;) is the error function. 



'27T 



r 2 + 7 



+ 8e" 



(A3) 



2. w = l/10 



P2(r) = 



e 5 



1000 r (r 2 + 45) + 5eio\/T07r 



25000tt (61tt + 48^) r 
+ 2875) + 20 (r 4 + 50r 2 + 175) erf (^j=) }• 



r r 4 + 190r 2 



(A4) 



3. cj = (35 - 2 V57) /712 « 0.0173 



In this case we obtain 



p 2 (r) = -u 3/2 p(Vujr), 



where 



c = 



7T 

384 



3 (67741 + 8855^) vr + 256 (l32 + 17y/&?) ^ (70 + 6y/57) vr 



(A5) 



(A6) 



and 



192 (-35 + 3V57) y 



178 



,_f8 (-34410443 + 4535711^) y 6 
35 -3^^ V J y 

+ 4 (-307865169 + 39250957^) y 4 + (l987335802 - 305824402^) y 2 

+ 1407537905^ - 1107674522l) + e^Vn(l6 (-34410443 + 4535711^) y f 

+ 32 (-394348481 + 50982845^) y 6 + 8 (-773680875 + 26538047^) y 4 

+ 56 (-358975185 + 8100493^) y 2 + 9 (-22356535675 + 2789329039^57^) 

+ Ae y2 * / — 178?r [l6 (-34410443 + 4535711^) y s + 32 (-85568903 
y 35 — 3y57 L V /V 

+ 10946667^) y 6 - 56 (-53901687 + 8766619^) y 4 + 72 (-272540317 

+ 33930025 Vdfjy 2 + 33 (-412010907 + 53643311^) er%) 1. (A7) 
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